########################################################################################
### Replication Code                                                                 ###
### Title: Overestimation of the Level of Democracy among Citizens in Nondemocracies ###
### Author: Eddy S. F. Yeung                                                         ###
### Version: March 5, 2022                                                           ###
########################################################################################

### Set-up ###
# Clean the R environment and set the working directory
# RStudio version 1.3.959 (macOS)
rm(list = ls())
setwd("~/Desktop/CPS_replication/figures") # save figures in a "figures" folder

# Load the required packages
library(haven)     # version 2.4.3
library(dplyr)     # version 1.0.7
library(sensemakr) # version 0.1.4
library(miceadds)  # version 3.11.6
library(extrafont) # version 0.17

# Import the dataset
df <- read_dta("~/Desktop/CPS_replication/dataset_final.dta")

# Rename the variables (for visualization purposes)
names(df)[names(df) == "growth_one_yr"] <- "GDP_Growth_Rate"
names(df)[names(df) == "new_rol"] <- "Rule_of_Law"
names(df)[names(df) == "new_gov"] <- "Govt_Effectiveness"

### Figure S19 ###
# Sensitivity analysis for Freedom of the Press
# Benchmark covariate: GDP growth rate
pdf(file = "Figure S19a.pdf", width = 6, height = 6)
model.temp1 <- 
  lm(diff_dem_vdem ~ new_fotp + university + female +
                     age + age_sq + married + unemployed + income +
                     social_class + ln_gdp + GDP_Growth_Rate +
                     Rule_of_Law + Govt_Effectiveness, 
     data = df)
sensitivity_1A <- 
  sensemakr(model.temp1,
            treatment = "new_fotp",
            benchmark_covariates = "GDP_Growth_Rate",
            kd = c(6, 12),
            ky = c(6, 12),
            q = 1,
            alpha = .05,
            reduce = T)
par(family = "Times", ps = 16)
plot(sensitivity_1A)
dev.off()

# Benchmark covariate: government effectiveness
pdf(file = "Figure S19b.pdf", width = 6, height = 6)
sensitivity_1B <- 
  sensemakr(model.temp1,
            treatment = "new_fotp",
            benchmark_covariates = "Govt_Effectiveness",
            kd = 1:3,
            ky = 1:3,
            q = 1,
            alpha = .05,
            reduce = T)
par(family = "Times", ps = 16)
plot(sensitivity_1B)
dev.off()

### Figure S20 ###
# Sensitivity analysis for Media System Freedom
# Benchmark covariate: GDP growth rate
pdf(file = "Figure S20a.pdf", width = 6, height = 6)
model.temp2 <- 
  lm(diff_dem_vdem ~ new_msf + university + female +
                     age + age_sq + married + unemployed + income +
                     social_class + ln_gdp + GDP_Growth_Rate +
                     Rule_of_Law + Govt_Effectiveness, 
     data = df)
sensitivity_2A <- 
  sensemakr(model.temp2,
            treatment = "new_msf",
            benchmark_covariates = "GDP_Growth_Rate",
            kd = 1:3,
            ky = 1:3,
            q = 1,
            alpha = .05,
            reduce = T)
par(family = "Times", ps = 16)
plot(sensitivity_2A, xlim = c(0, .61))
dev.off()

# Benchmark covariate: government effectiveness
pdf(file = "Figure S20b.pdf", width = 6, height = 6)
sensitivity_2B <- 
  sensemakr(model.temp2,
            data = df,
            treatment = "new_msf",
            benchmark_covariates = "Govt_Effectiveness",
            kd = c(5, 10),
            ky = c(5, 10),
            q = 1,
            alpha = .05,
            reduce = T)
par(family = "Times", ps = 16)
plot(sensitivity_2B)
dev.off()

### Figure S21 ###
# Sensitivity analysis for university education
# Benchmark covariate: income
pdf(file = "Figure S21a.pdf", width = 6, height = 6)
sensitivity_3A <- 
  sensemakr(model.temp1,
            treatment = "university",
            benchmark_covariates = "income",
            kd = c(2, 4, 6),
            ky = c(2, 4, 6),
            q = 1,
            alpha = .05,
            reduce = T)
par(family = "Times", ps = 16)
plot(sensitivity_3A)
dev.off()

pdf(file = "Figure S21b.pdf", width = 6, height = 6)
sensitivity_3B <- 
  sensemakr(model.temp1,
            treatment = "university",
            benchmark_covariates = "social_class",
            kd = c(5, 10),
            ky = c(5, 10),
            q = 1,
            alpha = .05,
            reduce = T)
par(family = "Times", ps = 16)
plot(sensitivity_3B)
dev.off()
